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Abstract 

Computing the atomic geometry of lattice defects — point defects, dislocations, crack tips, surfaces, or 
boundaries — requires an accurate coupling of the local strain field to the long-range elastic field. Periodic 
boundary conditions used by classical potentials or density-functional theory may not accurately reproduce 
the correct bulk response to an isolated defect; this is especially true for dislocations. Recently, flexible 
boundary conditions have been developed to produce the correct long-range strain field from a defect — 
efi'ectively "embedding" a defect in a finite cell with infinite bulk response, isolating it from either periodic 
images or free surfaces. Flexible boundary conditions require the calculation of the bulk response with the 
lattice Green function (LGF). While the LGF can be computed from the dynamical matrix, for supercell 
methods (periodic boundary conditions) it can only be calculated up to a maximum range. We illustrate 
how to accurately calculate the lattice Green function and estimate the error using a cutoff dynamical matrix 
combined with knowledge of the long-range behavior of the lattice Green function. The effective range of 
deviation of the lattice Green function from the long-range elastic behavior provides an important length 
scale in multiscale quasi-continuum and flexible boundary-condition calculations, and measures the error 
introduced with periodic boundary conditions. 

PACS numbers: 61.72.Bb, 61.72.Ji, 61.72.Lk, 61.72.Mm, 61.72.Nn, 62.20.-x 
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I. INTRODUCTION 



Lattice defects — e.g., interstitials, vacancies, dislocations, crack tips, free surfaces, interfaces, 
and boundaries — each play key roles in material properties,' and in order to understand defects, 
one must know their geometry. The far-field geometry for many defects is accurately described by 
anisotropic elasticity theory;^ However, the elastic solution often diverges near the atomic-scale 
center of the defect, and in many cases the center is difficult to investigate with current microscopy 
techniques. This is especially true of dislocations, which control plasticity in metals^ and can 
severely limit device utility in semiconductors.'' Only recently has the geometry and electronic- 
structure of an isolated dislocation been calculatedr^ this despite the rapid advances in computer 
hardware and density-functional theory methods. Previous density-functional theory calculations 
were limited by the long-range strain field of a dislocation which is incommensurate with peri- 
odic boundary conditions; hence, only dislocation dipoles*^-^ or quadrapoles'" " had been com- 
puted. The advent of "flexible" or "Green function" boundary conditions — first conceived by 
Sinclair et al.,^^ later redeveloped for crack propagation^^ and for dislocations and dislocation 
kinksii — made possible the relaxation of the core geometry of an isolated dislocation. For a re- 
view of density-functional theory methods applied to dislocations, see lisll . Flexible boundary 
conditions accurately treat the long-range strain field away from the defect by using the harmonic 
ideal lattice response in the form of the lattice Green function. The lattice Green function deter- 
mines the relaxed position of an atom given the force on it and its neighbors. Flexible boundary 
conditions have been used to model cracks,'^^*^^ dislocations and kinks in bcc metals with classical 
potentials cross-slip processes in fee metals,'^ isolated screw dislocations in bcc metals and 
ordered intermetallics with density-functional theory,i^i^ and even vacancies and free surfaces;— 
for a review of flexible boundary condition approaches to nanomechanics of defects, see jzil. 

Flexible boundary conditions are limited by the accuracy of the lattice Green function. Many 
closed-form results are known for the lattice Green functions of cubic lattices with nearest neigh- 
bor interactions.-^^ While the lattice Green function is intimately related to the elastic constants 
and dynamical matrix of a crystal, it has previously been computed for realistic potentials from 
relaxation of atom positions given an applied force. ^^ Rao et al. employed a "direct displacement" 
technique where separate relaxation calculations in a two-dimensional slab are used to numer- 
ically evaluate the lattice Green function for short range, while switching to the known long- 
range behavior of the elastic Green function.''' Woodward et al. used this same technique with 
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density-functional theory for Mo, Ta, and TiAl, and found the lattice Green function matched the 
long-range behavior at distances of only 5A, despite long-range metallic bonding.- However, this 
technique is dependent on the defect geometry — a lattice Green function computed for a [1 10]/ 2 
fee screw dislocation cannot be used for the fee edge dislocation with a threading direction of 
[1 12]/2. Moreover, relying on atomic relaxation can be prone to error in density-functional meth- 
ods when the applied forces become small. A more accurate and efficient approach instead relies 
on the dynamical matrix and elastic constants, which can be computed using standard techniques. 

What follows is a general and accurate method for the computation of the lattice Green function 
applicable for use in density-functional theory for a variety of defect geometries. In addition, we 
present and test an estimate of the error in the lattice Green function due to the geometry limita- 
tions of periodic -boundary conditions with density-functional theory. Currently available methods 
for computing the dynamical matrix in density-functional theory effectively produce a "folded" 
dynamical matrix, defined in an artificial supercell — whether they rely on an finite supercell or 
calculated on a discrete ^-point gri(i.^i21S^i2ZiS2i2^ However, the interactions in density-functional 
theory have an unknown range, likely to be larger than the artificial supercell. What is required to 
compute the lattice Green function is (1) a computational algorithm to accurately use the limited 
dynamical matrix information, and (2) an estimate of the error introduced from the dynamical 
matrix limitation. 

Section inireviews the harmonic response functions in a lattice — the dynamical matrix, and lat- 
tice Green function — and relation to continuum elasticity theory. Section Hill derives the general 
procedure for accurate numerical evaluation of the lattice Green function, with specific application 
for zero-, one-, and two-dimensional defects (point defects, dislocations, and boundaries, respec- 
tively). Section|lVlderives an error estimate for the lattice Green function using only the dynamical 
matrix computation from a single supercell and elastic constants. The error estimate is numerically 
tested using a simple-cubic lattice with random long-range interactions, and is shown to be accu- 
rate even with supercells far smaller than the interaction range. Finally, Section [V] concludes with 
discussion of applications to defect calculations and the inherent length-scales in quasi-continuum 
methods used in multiscale applications. 
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FIG. 1 : Schematic of harmonic lattice response. Displacement of atoms m in a crystal produce forces / on 
neighboring atoms; the forces are given by the dynamical matrix D for small displacements. Conversely, 
if an atom experiences a force, neighboring atoms must displace in order to accommodate the force; the 
displacement is given by the lattice Green function for small forces. 

n. HARMONIC LATTICE RESPONSE 

When atoms in a crystal are subject to applied or internal forces, they respond by displacing 
from their ideal lattice sites; and conversely, displacement from the ideal lattice sites produces 
internal forces on atoms. In the case of small displacements and forces, the atoms in the lattice 
respond harmonically. Harmonic response is characterized by a linear relationship between dis- 
placement and force, given by two different lattice functions: the dynamical matrix, and lattice 
Green function (Figure [T] shows the responses schematically). Below we review their definitions 
and connections to anisotropic elasticity theory. To simplify the notation, we assume a single atom 
Bravais lattice; however, the approach translates readily to multiple atom Bravais lattices. Ionic 
crystals have additional complexities that are not addressed here.— 

The dynamical matrix is well known from the classical and quantum theory of the harmonic 
crystal."^2i^ Let ^ and ^' be two lattice sites in a crystal, and i/(^) and it{^') the displacements of 



the atoms from their ideal lattice sites. Then we write the harmonic potential energy as 

M' 

where D{R - R') is a 3 x 3 matrix defined at lattice sites, and the double sum ranges over all lattice 
sites. For small displacements, the harmonic potential energy will equal the total potential energy 
^totai system up to a constant if the dynamical matrix components are 

^2f;total 



D AR - R') = 

dua{R)dub{R') 

In an infinite bulk lattice, there are several symmetry relations. First, _D is a function only 
of ^ - not ^ and ^' independently due to translational symmetry. Furthermore, D^,i,{^) = 
D^i^(-R) = Dj^^{R) due to inversion symmetry and independence of differentiation order. Finally, 
if we displace every atom identically, the total energy must remain a constant; hence, the sum rule 
^j}D{R) = 0. This sum rule has important consequences for the lattice Green function. 

The dynamical matrix linearly relates internal displacements and forces, and connects elastic 
strain to elastic stress through the elastic stiff"ness tensor. Given displacements u(M) at atom^, the 
internal forces /(^') produced at atom ^' are given by 

a T rharm 

f(R') = -—^ = -2]d{r- r')U(R). (1) 

An interesting special case of Eqn. ^ are displacements corresponding to a constant strain: u{^) = 
sR, where s is the strain tensor. The crystal response is a constant stress tensor a which is linearly 
proportional to the strain by Hooke's law: = Cs. This relationship is valid for small strains, and 
defines the fourth-rank elastic stiffness tensor C. Eqn. © gives a connection between the elastic 
constants Catcd and the dynamical matrix, 

- ^ D^il,(R)RcRd = V(Cacbd + Cadhc), (2) 
R 

where V is the volume of the unit cell. Eqn. © can also be derived using the method of long- 
waves.— Hence, the elastic constants contain information about long-range behavior of the dy- 
namical matrix. 

The lattice Green function**^ linearly relates the forces on each atom to its displacement from 
the ideal lattice. That is, 

il(M') = -J]GHR-R')f(R), (3) 



where ^(R-R') is the lattice Green function. It obeys similar symmetries to the dynamical matrix: 
G|^^(i?) = ^fj{-R) = ^a^^^- However, there is no sum rule for the lattice Green function.— 
At first, Eqn. © may not appear useful; to compute the forces in the harmonic potential, the 
displacements must already be known. However, if instead the forces on atoms in a simulation 
are computed using the total energy which is a function of relative atom positions, Eqn. Q 
allows one to relax the atoms to their ideal lattice positions. In that case, the displacements are 
not known when computing the forces. In particular, the lattice Green function is used to create 
flexible boundary conditions'^ where an isolated defect is surrounded by atoms that respond as 
if they are coupled to infinite bulk. This gives an accurate treatment of the long-range stress field 
of a defect (such as a dislocation) while using forces from close to the defect. 

As the long-range behavior of the dynamical matrix is connected to the elastic constants, 
the long-range behavior of the lattice Green function is connected to the elastic Green func- 
tion. The elastic Green function G is a continuum function that relates a force-field f{x) to 
the displacement-field u{y): 



The elastic Green function can be computed knowing only the elastic response of the continuum — 
i.e., the elastic constant tensor C^^{x) satisfies the partial differential equation 



where 6{x) is the Dirac delta-function. The lattice Green function must match the elastic Green 
function as 7? ^ oo, regardless of how long-ranged the dynamical matrix is. 

Lastly, the dynamical matrix and lattice Green function are inverses of each other. Substituting 
Eqn. ([T]) into Eqn. © gives 



where 6{R) is the Kronecker delta-function. Eqn. I© is not strictly solvable because D is singular 
due to the sum rule. The singularity is due to the lack of forces from a uniform displacement 
of all atoms; thus, the displacements from Eqn. © will be known only up to a constant shift in 
the entire lattice. This overall translational symmetry in the lattice Green function provides for 
the "flexibility" in flexible boundary conditions: bulk lattice response can be simulated without 
specifying an origin for the lattice. Mathematically, the singularity in D must be carefully treated 
to compute the lattice Green function accurately. 
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The computation of the lattice Green function is more tractable in reciprocal space. The lattice 
functions can be written as periodic functions of vectors k in the Brillouin zone (BZ) of reciprocal 
space,— 



In reciprocal space, the inverse equation Eqn. ^ simplifies to G^{k)D(k) = 1 for all k. The 

singularity of D is reduced to the gamma point k = 0, where D(0) = 0; for all other points, 

~, ^ — > -1 

G (k) = [D{k)] . The inverse is well-defined for metastable crystal structures; i.e., crystal struc- 
tures without unstable phonon modes. 

The computation of lattice Green function relies on accurate computation of the dynamical ma- 
trix. While computing the dynamical-matrix is straightforward for interactions with a finite cutoff", 
it is difficult for density-functional theory methods which may have long-range interactions (such 
as Friedel-oscillations). Two methods have emerged: direct force^^'^S'^^ and linear-responsci^^ 
Both methods compute the reciprocal-space dynamical matrix on a discrete grid of k-points in the 
BZ. This is equivalent to folding the real-space dynamical matrix into an artificial supercell. We 
use the folded dynamical matrix to compute the lattice Green function; thus, we need to evaluate 
the effect of the cutoff on the accuracy of the resulting lattice Green function. We do so with 
the elastic constants, which can be found separately by computing the response of a periodic cell 
to uniform strains. Eqn. © relates the elastic constants to the long-range behavior of the true 
dynamical matrix. This relation provides an estimate for the deviation of the long-range elastic 
Green function from the lattice Green function, which is turn gives an error estimate for using 
the folded dynamical matrix. More importantly, this estimate does not rely on a convergence test 
computation comparing increasingly larger supercells. 

Finally, it should be noted that the lattice Green function defined in Eqn. © can be modified 
for different bulk boundary conditions. Eqn. © defines in infinite bulk, called the 3D lattice 
Green function, and it is useful for computation of point defects. If the forces and displacements 
have periodicity along a lattice vector f, such as in a single straight dislocation defect, the 2D 
lattice Green function is used: Yjn G^(R + ni). Finally, if forces and displacements have periodicity 
along two lattice vectors and fj, such as in surfaces, grain boundaries and interfaces, the ID 
lattice Green function is used: Yjmn^(R + + "^2)- Despite the simple summations used to 
define the 2D and ID lattice Green functions from the 3D, the sums converge conditionally. It 
should be remembered that the "dimensionality" of the lattice Green functions refer to the degrees 




R 



BZ 
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of freedom for the lattice vector i? — remains a 3 x 3 matrix in all cases. The dimensionality of 
the defect (0, 1, or 2) plus the dimensionality of the lattice Green function (3, 2, or 1) sums to 3. 

III. COMPUTATION OF LATTICE GREEN FUNCTION 

The procedure for numerical computation of the lattice Green function separates the Fourier 
transform into pieces which can be inverse Fourier transformed accurately. The straightforward 
approach would be to discrete inverse Fourier transform the inverse of the dynamical matrix; 
however, this transform converges very slowly with increased grid spacing due to the second- 
order pole at the gamma point. The inversion of the dynamical matrix to compute the lattice 
Green function is still best performed in reciprocal space, where the large R behavior is exactly 
contained in the pole at k = 0. To accurately compute the lattice Green function requires an 
analytic treatment of the small k behavior separated from the rest of the Brillouin zone. 

The separation of the lattice Green function allows the inverse Fourier transform to converge 
by analytically treating the second-order pole. Moreover, the separation can be evaluated for any 

dynamical matrix, and for any dimension. The second-order pole in comes from the expansion 

— -> 

of D(k) for small k. 




R 




(6) 



The final expression is rewritten in terms of two functions of k of different order in k: k'^K^^\k) 
k'^A'-^Xk), where k = k/k. We relate the first function A^^^(^) to the elastic constants by Eqn. @, 




which gives A^^^(^) = y[^C^], where C is the fourth-rank elastic stiffness tensor. On the other 
hand, the quartic function A'^'^^(^) has no similar simple connection. With the definitions of A^^^ 
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and A^'^^ the lattice Green function expands for small k as 



G^{k) = {Dm' 



= [k^A^^\k) - k'^A^^k) + 0(k^)] 



-1 



= r^[A(2)(^)]"'[l - k^A^'^\k)[A^^\k)]'^ + Oik"^)] ' 
= r2[A(2)(^)]"' + [A^^\k)]'^A^'^\k)[A^^\k)]~^ + 0{k^) 



= G\k) + G^{k) + 0{k'), 



dc/ 



where 



-1 



G^{k) = k-\A^^\k)\ 

1 ^ ^ -I 
[kCk-] , 



(7) 



and 



C^^k) = [A^^\k)] 'A^^\k)[A^^\k)]^ 



i 



G^(k). 



(8) 



The first function is the second-order pole at the gamma point, which is the Fourier transform 
of the elastic Green function. The second function G'*'^ is independent of |^|, representing a discon- 
tinuity at the gamma point in the lattice Green function that appears only after the second-order 
pole is subtracted out; this function is called the discontinuity correction. 

This expansion is used to separate the Fourier transform of the lattice Green function in the en- 
tire Brillouin zone into three pieces: the elastic Green function, discontinuity correction, and semi- 
continuum correction. We introduce the continuous and difi'erentiable cutoff function /cut(^/^max) 
with parameter < or < 1, 



.^ut(-^) 



1 



3(f^)" - 2( 




<x <a 
or < X < 1 > 

1 < JC 



(9) 



where fcmax is the radius of a sphere inscribed in the Brillouin zone. While final evaluation of ^(^) 
is independent of a, all computations to follow use a = 1/2. Then, the semicontinuum correction 
is defined for k in the first Brillouin zone as 



~ -1 



G''(k) = [D(k)] ' - (G^(k) + G^\ky)fUklk^^) 



(10) 




lattice GF elastic GF discontinuity semicontinuum 

correction correction 



FIG. 2: Separation of lattice Green function for a square lattice in 2-dimensional reciprocal space into elastic 
Green function, discontinuity correction, and semicontinuum correction (note different vertical scales). The 
lattice Green function has the periodicity of the reciprocal lattice, and a second-order pole at the gamma 
point. The elastic Green function scales as k~^, and is cutoff to smoothly go to zero at the Brillouin zone 
edges. The removal of the second-order pole creates a discontinuity independent of \k\ at the gamma point; 
the discontinuity correction removes the discontinuity and smoothly goes to zero at the Brillouin zone edges. 
The remaining difference between the lattice Green function and the first two terms is the semicontinuum 
correction, which is smooth everywhere in the Brillouin zone. 

The elimination of the second-order pole at the gamma-point by using a cutoff version of the elastic 
Green function is related to the semicontinuum method of Tewary.-^ However, his semicontinuum 
approach used a Gaussian cutoff which does not vanish at the Brillouin zone edge, and does not 
treat the discontinuity produced at the gamma point. The final lattice Green function is the sum of 
three pieces: elastic Green function, discontinuity correction, and semicontinuum correction. 

Figure|2|shows an example of the separation of the lattice Green function into the three terms for 
a square lattice. The lattice Green function shown comes from a square lattice with lattice constant 
qq = n and nearest-neighbor interactions. For this case, G^(kx, ky) = (sin^(;r/r^/2) -i- sin^(;r^y/2))"^ 
The second-order pole at origin is given by the elastic Green function G^{kx,ky) = Al{n\k\f; it 
is multiplied by the cutoff function with fcmax = 1 so as to vanish at the Brillouin zone edge. 
Subtracting the pole from the lattice Green function produces a function with a discontinuity at the 
gamma point. The discontinuity at the origin is given by the discontinuity correction G^^{kx, ky) = 
(k^ + k^,)/{3\kf) which is multiplied by the cutoff function. Subtracting the discontinuity produces 
the semicontinuum correction, G^^(k^, ky), given by Eqn. (flUI) . 

The evaluation of the lattice Green function in real space is accomplished by inverse Fourier 
transforming the semicontinuum correction G^^, the cutoff elastic Green function G^/cut, and the 
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TABLE I: Overview of lattice Green function computation for different dimensionality. The dimensionality 
of the lattice Green function is determined by the type of defect being simulated: the defect dimensionality 
plus the lattice Green function dimensionality is three. While the lattice Green function has the same form in 
reciprocal space, the periodicity determines the range of Brillouin zone integration, and the functions used 
to expand the k dependence of the elastic Green function and discontinuity correction. The range of BZ 
integration produces different large R behavior for both the elastic Green function — also given by elasticity 
theory — and the discontinuity correction. The one-dimensional case has no k dependence, so there is no 
angular expansion nor is a discontinuity correction required. 



Defect type: 
(dimensionality) 


3D 

point (OD) 


2D 

dislocation, 
crack tip (ID) 


ID 

free surface, 
boundary (2D) 






plane(s) ± to 




Brillouin zone 




line(s) ± to 


full BZ 


threading 


integration: 






surface plane 




direction 


Angular 


spherical 






Fourier series in 






harmonics 




N/A 


expansion in 




plane e'"* 


Brillouin zone: 


yim(Sk,<l>k) 








Large R elastic 




-lnR + R° 


R 


Green function: 








Large R 

discontinuity 

correction: 


R-3 


R-2 


N/A 



cutoff discontinuity correction G f cut- The semicontinuum correction G^^(k) is evaluated on a 
discrete k-point grid in the Brillouin zone using Eqn. dTOt : inversion of the dynamical matrix for 
small k must be handled carefully to avoid numerical noise. A discrete inverse Fourier transform 
converges well with grid spacing because G^^ is smooth throughout the Brillouin zone. The cutoff 
elastic Green function G^/cut and discontinuity correction G'^^fcut are expanded as functions of k 
using spherical harmonics or a Fourier series depending on the dimensionality of the problem. 
In this form, their inverse Fourier transforms can be analytically reduced to a one-dimensional 
integral of non-singular functions over a finite range that is computed numerically to the desired 
accuracy. The details of this reduction depends on the periodicity of the lattice Green function. 
TablelUgives a brief overview of the results, and Table HH gives a summary of the equations. 
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A. 3D lattice Green function: OD defects 



To facilitate inverse Fourier transformation, the elastic Green function (Eqn. dTJl) and discon- 
tinuity correction C^^ (Eqn. ©) are expanded as a spherical harmonic series, whose coefficients 
are computed numerically. The expansions 

1=0 m=-l 1=0 m=-l 

are truncated for / > L^ax; ^max is chosen for each expansion so that the Im components above 
Lmax are less than 10"'^ of the largest Im component below L^ax- Moreover, as both and G^'^ 
are symmetric with respect to inversion, only even I values are nonzero. The normalized spherical 
harmonics are given by 



V An (/ + m) ! 

for k = (sin cos 0, sin 6 sin 0, cos 6), where P'J'ix) is the associated Legendre polynomial without 
the (-1)™ phase.— To compute the spherical harmonic expansion, the elastic Green function and 
discontinuity correction are evaluated on a spherical N x N grid \k\ = 1 given by 0, = Ini/N, and 
6j = arccos(My), where uj are the roots of the A^"^ order Legendre polynomial Pn{u). This grid 
allows the computation of expansion elements up to L^ax = N/2 - 1. The spherical harmonic 
components are evaluated by (1) discrete Fourier transforming the 0, grid to m components, then 
(2) using Gaussian-Legendre quadrature with weights Wj on the Oj grid to produce / components: 



- 2^ V p/+l(/-m)! Y ^~i""P.i=^(n ^ ^ 

and identically for G'^'^. As a final step, it is useful to reduce numerical error in the expansion by 
explicitly enforcing the point group symmetry of the lattice on the expansion; this is done using 
the Wigner D-matrix,^^ modified to take into account the effect of a symmetry operation on the 
3x3 matrix elements. 

Given the spherical harmonic series, inverse Fourier transformation reduces to a single integral 
over a finite range. Writing the inverse Fourier transform integral in spherical harmonics over the 
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BZ gives 

GH^) = Gf,n-^ j dk Uik/k^,.) Jj d^k e''''^'-%^(k), 

''" Att 

J ^max 

^\^) = J^G^^— J dk k^U(k/k^^) J J d^k e'"''^'-^^YUk), 



by virtue of the cutoff function. The double integral over k is evaluated analytically as Eqn. (IA3I) . 

so 

T ^max 

G^(M) = Gf,„Y,JR)(-lf^^ j dk fUklk^,,)ii{kR\ 

Im ^ Q 

T ^max 
^max T 7- /-» 

G^\M) = Y GtYimm-lf^^_ J dk k^Uik/k^,,)j,m, 



Im Q 



where ji(x) is the spherical Bessel function. 

The finite one-dimensional integrals are smooth functions that can be evaluated numerically to 
required accuracy. For the special case of 7? = 0, 

I. 

Jdk fcut(k/km^)ji(0) = Sifikn 



^max 

I i^fc /cut(^/^max)j/(0) = 5/,ofc„ 



30 



For R i^O,v/e define ff'Xx) and fi^\x) as 



2 

//°\.«) = - \ du ji{u)f^ut{u/x), 
^ Jo 



r(2)/ 



//^^x) = - f du u^ji{u)fcutiu/x), 

n Jo 



(11) 



so 



^ J^fc U{klk^,,)h{kR) = l^f\k^,,R), 





i. 

"-max 



^ J dkk^fUk/k^MkR) = ^f\k^^R). 



13 



The // functions are evaluated numerically by splitting the integrals into intervals between roots of 
ji{x), and then using the QAG adaptive integration algorithm with 61 point Gauss-Kronrod rules 
from QUADPACK.— An important limiting case is for R ^ oo where the functions can be evaluated 
analytically. From Gradshteyn and Ryzhil^ expression 6.561.14 gives for even / 



These results are summarized in Table HH 

The inverse Fourier transform of the semicontinuum correction G^^ is performed with a discrete 
transform on a grid in the Brillouin zone. There are different techniques for constructing a k-point 
meshr^i^ but a uniform grid of ^-points centered at the gamma point inside the BZ suffices. The 
primary requirement is that each k-point lie in the first BZ; and in that way, the points given by 
1^1 < ^max form a sphere. The spacing of the grid is determined by the largest magnitude lattice 
vector 7?max in the desired domain of G^(^). To avoid aliasing errors, the grid spacing Ak must be 
smaller than In/R^i^, though a smaller spacing is preferable. For large R, substituting the elastic 
Green function for the lattice Green function introduces only small errors, hence reducing the 
effective R^ax and k-point mesh that is used. We estimate the deviation in detail in Section Hvl 

B. 2D lattice Green function: ID defects 

The introduction of a threading direction reduces the lattice Green function to a two- 
dimensional slab and modifies the inverse Fourier transformations. The forces and displacements 
of atoms around a dislocation line or a crack tip have a periodicity given by a threading lattice vec- 
tor L The periodicity is represented in the lattice Green function by the 2D lattice Green function, 
2n + n?). As with the 3D lattice Green function, evaluation of the 2D lattice Green function 
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is best performed in Fourier space, and inverse Fourier transforming to real space. Then, 

CO 

n=-oo 

= t^JJ!^'''-"'-'"'''^^^^ (12) 



BZ 



^lieBZ BZ 

where the (finite) summation is over k\\ = 2nmtl\t\^ (m integer) that are inside the BZ, and two- 
dimensional integration is over kj^ that are perpendicular to t and inside the BZ. This is by virtue 
of the summation over n which produces a delta function on exp(/fc ■ i) - I. Eqn. dT^ still has a 
pole in G to contend with, but it lies purely in the plane of k\\ = 0. Hence, for k^ 0, the value 
of G^ = [D] ^ is used, and a discrete inverse Fourier transform is performed. Then, the remaining 
difficulty is the l/k\ pole at the gamma point in the 2D inverse Fourier transform. 

The pole at the gamma point in 2D is split into three contributions for inverse Fourier transfor- 
mation: elastic Green function, discontinuity correction, semicontinuum correction. The elastic 
Green function G^ and discontinuity correction G'^^ are expanded as a truncated Fourier series in 
the plane of kj^, 

-L n=0 n=0 

where cpt is the angle of kj^ relative to an (arbitrary) normalized in-plane reference direction «x 
(hj_ ■ t = 0). The truncation A^^ax is chosen for each expansion so that the n components above 
A^max are less than 10"'^ of the largest n component below A^max- Since both G^ and G'^^ have 
inversion symmetry, only even n values are nonzero. The Fourier series components are evaluated 
by computing G (k) onaN element circular grid l^^l = 1 at a series of angles 0,- = Ini/N relative 
to Hx- The discrete Fourier transform gives 

. N-l 

G^n = -J]e-'^'''G^{hm, 

i=0 

and identically for G'^'^. Note that G^{k) and G'^^ik) are the same functions that appear in the 3D 
lattice Green function (given by Eqn. Q and Eqn. dUl); for the 2D lattice Green function, they are 
only evaluated in the plane through the gamma point. 

Given the Fourier series, inverse Fourier transformation reduces to a single integral over a finite 
range. The terms of Eqn. dT^ have no singularities, so they can be evaluated numerically 
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using a discrete inverse Fourier transform with a discrete grid for k^^, where construction of this 
grid is described below. Hence, the elastic Green function and discontinuity corrections are only 
evaluated for k\\ = 0. The inverse Fourier transform integral over the BZ gives 



G\R) = 0] = 

_ cos(0i— 



'- " max T7 r* y^i 

^ \t \{lnY J k J 



where Rj^ = sJr'^ - ■ is the magnitude of ^ perpendicular to t, and (pa is the in-plane angle 
of R relative to hj_ {<pR = arccos((nx ■ R)/R±))- The integral over (f>k is given by expression 8.41 1.1 
in Gradshteyn and Ryzhik— as 2n(-l)''^^Jn{kRj_) expiin^R) where Jn{x) is the Bessel function, so 

\T ^max 
^ Vmax T 7 /-* 

GE(^) = y G^^^"*(-l)"/2 dkk-'U(k/k^M(kRJ, 

If 

\j '*-max 
^jjTiax T 7 

G'\R)=yGf^e'"'^'^i-lf'-— dkkU(k/kmM(kRJ. 

The finite one-dimensional integrals are smooth functions for n > that can be evaluated 
numerically to required accuracy. For the special case of i?x = and n 0, the integrals over k 
are zero. For n ^ and Rj^ 4^ 0, we define i^J,°'(jc) and i^i'^'(x) as 



so 



F^^\x) = I duU ^Jniu)fcutiu/X), 

Jo 

Pn^{x) = I du uJn(u)fcut(u/x). 

Jo 



dk k-'fUk/k^,MkRJ = ^Fr(^maxi?x), 





(13) 



V r _ ^ (2) 

- I dk kfcut{k/kyj^i^^)J„(kRj^ = - -rrF^ (kmaxR±)- 

In J InRj^ 



As in the 3D case, the F„ functions forni^O are evaluated numerically by splitting the integrals 
into intervals between roots of Jn(x), and then using the QAG adaptive integration algorithm with 
61 point Gauss-Kronrod rules from quadpack.— For u < 10"^ in the integral, using J„{u)/u « 
l/2(u/2y^^ eliminates the division by zero. An important limiting case is for Rj^ ^ oo where 
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F,'\k R) 

max ^ 




7? 



FIG. 3: Evaluation of the R scaling for the circularly symmetric portion of the 2D elastic Green function 
given by i^|)'\^max^)> with a = 1/2 and /c^ax = ^- The long range behavior of the elastic Green function 
in two dimensions scales as - ln{R). The cutoff function retains the correct large R behavior, with small 
deviations after R = 2m lattice units. Moreover, the cutoff removes the divergence at the origin. 



the functions can be evaluated analytically. Expressions 6.561.14 and 6.621.4 in Gradshteyn and 
Ryzhik''^ give for n ^ 



1 



lim F^;;'>(x) = -, lim F)f\x) = n. 

The finite one-dimensional integrals for the n = case require additional analytic manipula- 
tion to be evaluated numerically. Figure |3] shows the convergence to the long-range behavior for 
F''o\kmiixR±)- The function Ff\x) is well-behaved for all values of x, and the limiting case of 
X ^ oo is as given above. The Rj^ = 0, n = integral for the discontinuity correction is 



dk /:/cut(^/^max)-/o(0) = kl 



1 (1 - aX3a + 7) 

2 20 



The n = integral for the elastic Green function does not converge as written, because the limit 
as i?x ^ diverges. To perform the integration, it is useful to remember from elasticity theory 
that the 2D elastic Green function in real space scales as ln(7?x)- The Fourier transform of In |Pl in 
two-dimensions is 



d'r e'^^-^lnl^ 



dr r\n(r)Jo(kr) = 
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which means the inverse Fourier transform integral is 

CO 



Using this relation, the pole at A: = can be evaluated analytically as 

"■max 

Cdk 





j Y UkRJ + Jy ■^"^'^^^^ [/cut(^/^max) - 1] 



^, / afemax \ 1 y (-1)" / gfemax-^X 

" 2 ; ^ 2Zj(„ + l)(n + l)!2^ 2 I 

'^max 

+ I Y Jo(kR±)fcut(k/k^^^), 

where y is the Euler constant (y 0.5772156649). The remaining integral is evaluated numeri- 
cally as before, and the series is summed numerically to within 10"' ^ The expression is finite for 
R± = 0, and in the limit i?x ^ oo becomes - ln(i?x) (cf. Figure|3l). The R± = 0, n = integral for 
the elastic Green function is 

^max 

Cdk 

I Y /cut(^/^max)^0(0) = In(^niax) + - In 2] 



6a\3 - a)lnQ' - (1 - a)(5a^ - 22a + 5) 
^ 6(1 -q;)3 ■ 

These results are summarized in Table HH 

The inverse Fourier transform of the semicontinuum correction G^^ is performed with a discrete 

transform on a grid lying on planes in the Brillouin zone. The planes are specified by the threading 

direction in the lattice t, to form a planar grid requires two in-plane lattice vectors and m^. 

All three vectors are mutually perpendicular, though not normalized. The N x M grid is the 

combination of k\\ and kj^, with 

J* 27rf 27Tnj_ n 2nm^ m 
kit, n, m) = -—t + ^ ^ 1 — - — , 
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where t, n, m are integers ranging over the interior of the BZ. The integers and M specify the 
in-plane grid spacing, and must be chosen sufficiently large to remove aliasing effects out to i?max- 
As for the 3D lattice Green function, the deviation between the 2D lattice Green function and 2D 
elastic Green function decreases with distance, thus requiring the computation of the lattice Green 
function out to a fixed distance dependent on required accuracy. 



C. ID lattice Green function: 2D defects 



The introduction of an infinite surface or boundary reduces the lattice Green function to a 
one-dimensional column and modifies the inverse Fourier transformations. The forces and dis- 
placements of atoms away from a boundary — be it a free surface, grain boundary, or interface — 
have a periodicity given by two non-parallel lattice vectors fi and ?2 lying in the boundary plane. 
The periodicity is represented in the lattice Green function by the ID lattice Green function, 
Yjinn^i^ + f^h + nti). As with the 3D and 2D lattice Green functions, evaluation of the ID 
lattice Green function is best performed in Fourier space, and inverse Fourier transforming to real 
space. Then, 

oo 



ni,n2=~oo 



BZ 



KplaneSBZ BZ 

where the (finite) summation is over 



^plane — 27r- 



(mi and m2 integer) that are inside the BZ, and one-dimensional integration is over kj^ that are 
parallel to fi x ^2 and inside the BZ. This is by virtue of the summation over ni and nj, similar to 
the 2D lattice Green function. Eqn. still has a pole in to contend with, but it lies purely 
on the line where A;piane = 0. Hence, for fcpiane ^ 0, the value of = [D] is used, and a discrete 
inverse Fourier transform is performed. Then, the remaining difficulty is the l/k^ pole at the 
gamma point in the ID inverse Fourier transform. 
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The pole at the gamma point in ID can be spUt into two contributions for inverse Fourier 
transformation: elastic Green function and the semicontinuum correction. For one-dimensional 
variation along ^j^, the elastic Green function is 

where the factor = A*^^^ (t[ xpi/\h x hlj depends on x and there is no remaining dis- 
continuity at the gamma point. Thus, the semicontinuum correction no longer vanishes at the 
gamma point, but instead smoothly approaches a constant value. Thus, the only piece to be treated 
analytically is the l/kj_ pole at the origin. 

The inverse Fourier transformation of the elastic Green function requires the evaluation of a 
single integral. The elastic Green function in real space is 

GE(^) [^11 = 0] = — ^ r ^ e-''^--^G^(k^)U(k/k,^) 

- t, X to J 271 



\h X t2\ 

BZ 

k, 

\h X t2\ 

^max 



J 2n' 



where is the (positive) magnitude of R perpendicular to the plane given by t\ and ?2- As with 
the 2D elastic Green function, the pole is separated off and related to the known Fourier transform 
giving 

^max ^max 

J* 2ji ^ fcutik/kraax) ~ J^^^ /cut(^/^max) 

^max ^max 

f cosjkRJ r cosiJcRj^ 

^ I ni? I ni? (/cut(K//Cmax) - 1) 



ak„ 
1 





/COS(^i?x) 
dk (/cut(^/^ax) - 1) 

l^^l CV7 J3 ^ cos(/:i„ax^±) 

>3l(,Kmax«±J 7 

^max 

r cosjkR^) 

I (/cutiK/fcmax) - 1). 



max 
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where Si(x) = sm{t)/t dt is the Sine integral. The remaining integral can also be evaluated in 
closed form, but the expression is lengthy. Two important values are i?x = and R± oo, which 
are 

hm G^(R) = - — —— — - 

lim G^{R) = --\RA 



The limiting behavior of G^ ~ \Rj_\ from elasticity theory is recovered. These results are summa- 
rized in Table Inl 

The inverse Fourier transform of the semicontinuum correction G^'^ is performed with a discrete 
transform on a grid in lines through the Brillouin zone. The grid spacing along the line must be 
sufficiently small to remove aliasing effects. As with the 3D and 2D lattice Green functions, the 
deviation between the ID lattice and elastic Green functions decreases with distance. Thus, the 
elastic Green function may be substituted at a fixed distance, and requiring the computation of the 
full lattice Green function for a finite set of points. 



IV. ERROR ESTIMATION FOR LGF 

Table HIl shows that as R becomes larger, the lattice Green function asymptotically matches 
the elastic Green function; this matching provides the basis for an error estimate of the lattice 
Green function. The elastic Green function can be computed knowing only the elastic constants; 
in turn, the elastic constants can be computed even for interactions without a fixed cutojf, such 
as density-functional theory. Hence, while the dynamical matrix computational may induce an 
artificial cutoff, the asymptotic limit of the lattice Green function is known exactly. Then an 
estimate of the error in the lattice Green function can be determined by estimating the deviation 
between the elastic Green function and lattice Green function. Surprisingly, an accurate estimate 
can be obtained using the elastic constants and the dynamical matrix from an artificially folded 
supercell, even if the dynamical matrix has non-zero elements outside the supercell. Hence, a 
single approximate computation of the dynamical matrix in a supercell together with the elastic 
constants provides an estimate of the accuracy of the supercell computation. This is shown in 
detail below. 



21 



TABLE II: Summary of equations for lattice Green function computation for different dimensionality. The 
split of the lattice Green function into three pieces is given for each, along with the angular expansion. The 
lattice Green function in real space, the limit of large R, and value at 7? = are also given. The cutoff 
function /cut and parameter a are defined in Eqn. Q; ^max is the radius of a sphere inscribed in the Brillouin 
zone (BZ). All k are restricted to be inside the first BZ. The finite BZ summations are done over a grid 
of A'^kpt points. The function A(x) is 1 for ;c = 0, and elsewhere. The 3D integrals ff\x), ff'^ix) and 
the 2D integrals Ff\x), F^n\x) are defined in Eqn. (fTTl and Eqn. (I13t . respectively. For the 2D case, the 
periodicity is defined by a threading lattice vector t, and /?x is the magnitude of R perpendicular to t. In the 
ID case, the periodicity is defined by two non-parallel lattice vectors ti and ^± is the magnitude of R 
perpendicular to the plane of Ti and Pi- 
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FIG. 4: Connection of solution of continuum differential equation mapped onto a lattice equation. The 
continuum differential equation that defines the solution on the left can be discretized by introducing a grid, 
and approximating derivatives with finite differences over the grid to produce a lattice equation. In the limit 
that the grid spacing becomes small compared to the length scale of variation of the solution, the discrete 
approximation matches the continuum solution. This mapping can also be reversed, by starting with a grid 
and a lattice equation, then taking the limit of zero grid spacing to produce a corresponding continuum 
differential equation. 

A. Derivation 

The asymptotic connection between the lattice Green function and the elastic Green function 
can be understood by viewing the lattice Green function as a "numerical grid" solution to the elas- 
tic Green function differential equation, as in Figure IH The mapping of a continuum differential 
equation onto a discrete grid with a lattice equation is a well known method for the numerical 
solution of multidimensional partial differential equations. The (partial) derivatives can be ap- 
proximated using finite differences on the grid. As the grid spacing becomes small compared to 
the length scale of variation of the solution, the continuum solution is recovered. Moreover, this 
mapping can be reversed: given a lattice equation, taking the limit of zero grid spacing can recover 
the continuum partial differential equation. In the case of the lattice Green function, the grid is 
defined by the crystalline lattice, the lattice equation by Eqn. Q and corresponding continuum 
differential equation by Eqn. ©. 

The analogy of the numerical solution of partial differential equations provides the basic idea 
for the estimation of the deviation between the lattice Green function and elastic Green func- 
tion. In finite difference applications, an estimate of the discretization error can be determined by 
substituting the true continuum solution into the discrete equation, and using Taylor series to ap- 
proximate the deviation;^ For the lattice Green function equation, it is the elastic Green function 
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that is an approximation, but the methodology for error estimation is identical and provides the 
deviation between the lattice and elastic Green functions. 

The relative deviation of the lattice Green function from the elastic Green function can be 
extracted using the real-space lattice Green function equation Eqn. Q. We begin by defining the 
relative deviation e'-'^(^), for 7? > 0, as 



Note that while G^(^) is only defined at lattice sites, G^(^) and e'^^(^) are continuum functions. 
To simplify this expression, we introduce the zeroth and second order deviation functionals A^°^[] 
and A'^^^[] of a continuum function f{R) around a point R, 



These functionals describes the deviation between the second-order finite diff"erence expansion of a 
function and the true value, with or without using the Taylor expansion. For small x, A'^*'^[/](^, x) « 



^x^ df/dR^ and A^^\f](^, x) ~ ^x"^ df/dR\ Since = D{-x), A'^\] and A'^\] can be used 




and substituting into Eqn. Q for 7? > to get 




A(2)[/](^, f) =l[f(R + f) + /(i - f)] 
-lf-VV/(^)-f. 



f(R) 



to writes^ 





[gH^)s^^(R) + A^^\G^e^^]{R, f)] 



Substituting these into the lattice Green function equation gives 
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FIG. 5: Separation of error estimation summation Eqn. (fTSl into regions 1 and 2. The different regions can 
be simplified by using different approximations for A*^"^ [G^e*-"^] and A^^^[G^]. In region 1, jc » /?, and the 
deviations scale as -G^{R)e^^{R) and /R^G^{R), respectively. For region 2, x «; /?, and the deviations 
scale as /R^G\R)e^^{R) and x'^/R'^G\R), respectively. 

The first term is zero because of the dynamical matrix sum rule. The second term is simplified to 
zero by using Eqn. Q, 

X abc 

= J](C,ajh + C,tja)^a^t^^{^) = 0, 

abc 

which is zero by virtue of Eqn. ^ and applying the interchangeability of partial derivatives, sym- 
metries of the elastic Green function and elastic constants: = G^., Ccajb = Cjbca = Cjbac, 
Ccbja = Cjacb = Cjabc, and VflVfo = VfoV,,. Thus, 

^ A^2)[G^](^, f)D(f) + A^°>[G^e°^](^, f)D(f) = 0. (15) 

X X 

To evaluate e^^(^) in Eqn. (fT51) requires approximations for hP\^^^^ and A^^^[G^]. These 
are built by using the x R and x :s> R asymptotic values over the two regions in Figure |5l 
For region 1 (x » R), both G^e°^(^ + x) and G^s°^(^ - f) are much less than G^{^)s^^{^) 
in A^°\^s°^]{R,x). For region 2 (x ^ R), we use the fact that G^(7?) ~ l/R and assume that 
^^{R) ~ 1 for some power /3 to be determined. Using the Taylor series for small x and 
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choosing the maximum x ■ R = I gives the approximation in regions 1 and 2 



-1 .x>R 



Similarly, for region 1 (x » R), A*^^^[G^] is dominated by the quadratic growth of the last term; as 
G^(i?) ~ 1 /R, this gives x^ /R^^{R). Using the next order in the Taylor series for small x in region 
2 «; i?), with the maximum value of x ■ R = 1 gives the approximation for regions 1 and 2 



\x^/R^ .x>R 
[x^/R^ ■.x<R. 
Substituting these approximations into Eqn. (fTSt gives 



[ |Al<fi li\>R 

= G-(R)Uj]x^m^^,J]^m\- 

[ \-^<R |.v1>ff j 

The two sums over \x\ > R can be simplified by using (1) the sum rule for the dynamical matrix, 

-_^D(f) = +_^D(f), 

|-fl>-R |i1<R 

and (2) using the elastic constants with Eqn. @, 

.v>|^| ^ f<|^| 

^^^Caxhx ~^ Cayby ^azbz 
- Caxhx(R) - Cayby(R) - Cazhz(R)\ 

^ ^V6C^,iR), 

where Cahcd(R) refers to the elastic constants derived using Eqn. ©, summing only over lattice 
sites X < R. We will show shortly that £^^(R) ~ R~^, so after dividing out G^(^) and substituting 
yS = 2, we have 



[ |i1<R k1<R j 

1 4 ^ 2V 



(16) 
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Eqn. (O can be further simplified from the small and large R limits. In the small R case (large 
region 1), the first terms of the left-hand and right-hand sides of Eqn. (fT6b become negligible. For 
small R, the dominant piece of l]|f<R D(x) is D(0), so 

In the large R case (large region 2), the second terms of the left-hand and right-hand side of 
Eqn. (O become negligible. The summations Jji^<rX^D(x) and Tj\j^<R x'^ D{x) are related to the 
/ = spherical harmonic expansion of and G'^'^ (c.f. Appendix |B] and Eqn. (|Bl|l '). so 

rpainn 9 ^ -' 



3R^ 

^00 



where C^^iR) is evaluated using the truncated dynamical matrix. Note that both pieces scale as 
R~^, justifying the earlier choice of jS = 2. Because the region 1 estimate falls off faster as R 
becomes large, and the region 2 estimate goes to zero as R goes to zero, the final approximation is 
to sum the two pieces 

rp 2V , 10 &R) 

The main feature of Eqn. (fTTl) is that the two region estimates can be determined using a single 
supercell calculation, even if the dynamical matrix lacks a finite interaction cutoff. The region 1 
estimate is computed by comparing the true elastic constants to the elastic constants computed 
from a folded dynamical matrix, and using the supercell dimension for R. The dynamical matrix 
can result from a direct force computation in a finite supercell, or evaluating D for a finite k-point 
grid. In the latter case, the inverse grid spacing provides the value for R. The region 2 estimate is 
computable because G^^ is only summed over x<Rfor D_{x). In both of these cases, it is assumed 
that the effect of folding the dynamical matrix into the supercell is approximately equivalent to 
truncating it outside the supercell. 



B. Numerical example of error estimate: simple-cubic lattice 

While Eqn. dTTl has the advantage of being computable for long-ranged dynamical matrices, 
it is not clear if too much accuracy has been lost in the series of approximations; so a numerical 
example is used to highlight the range of applicability. A series of pseudo-random long-range 
dynamical matrices are generated on a simple-cubic lattice with characteristics related to real 
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material systems; and for each, the lattice Green function, relative deviation to the elastic Green 
function, and region 1 and 2 estimates are computed. 

The dynamical matrices are generated using D(R) ~ sm(nR/ao)R~'^, cutoff at 25ao, with lattice 
constant cio = 1. The functional form is chosen to provide a long-range interaction, whose falloff 
is still fast enough to produce finite elastic constants in Eqn. ©. The sin(;7ri?/ao) functional form 
produces a Friedel-like oscillation, as might be expected in a metallic system. The dynamical 
matrix elements at each site are pseudo-random numbers from a Gaussian distribution with mean 
and standard deviation sm(7rR/ao)R~'^. The dynamical matrix is symmetrized using the cubic 
point group. The elastic constants and phonons are computed; if there are unstable phonons, or 
the elastic anisotropy is greater than 3, the dynamical matrix is rejected. 100 random, stable, 
long-range simple-cubic dynamical matrices are generated in this manner; for each, the lattice and 
elastic Green functions along with the relative deviation are generated. The dynamical matrix is 
"folded down" into supercells from 2 x 2 x 2 to 20 x 20 x 20 to compute the region 1 and 2 estimates 
in Eqn. (fTTt . 

Figure |6l shows the true deviation and estimates from our test case for both a single example, 
and the average results from the 100 dynamical matrices. As expected from the derivation, the 
region 1 estimate dominates for small R, and falls off as the supercell becomes large enough to 
accurately produce the elastic constants. The region 2 estimate becomes important for large R, 
capturing the long-range effect from the discontinuity correction. What is especially encouraging 
is that the error estimate is accurate even for small supercells — such as 2 x 2 x 2 — where the 
supercell dynamical matrix calculation is clearly inaccurate due to the long range. This is perhaps 
the most impressive feature of Eqn. (flTl : Even when the dynamical matrix calculation comes from 
a small supercell, the known elastic constants can still provide an accurate error estimate without 
requiring comparisons to larger supercells. Hence, a supercell- size effect estimate on the lattice 
Green function computation is provided from a single supercell dynamical matrix computation. 

V. DISCUSSION 

The deviation between the lattice Green function and elastic Green function in Eqn. (fTTl can be 
described by a single length scale /eias that characterizes the recovery of continuum elastic behavior 
from atomistic lattice behavior: £^^(R) ~ (/eias/^)^- This length scale determines the range out to 
which the lattice Green function should be computed in lieu of the elastic Green function. For 
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FIG. 6: Relative deviation between EGF and LGF for simple-cubic test case. The points give the deviation 
between the lattice Green function computed with the full dynamical matrix and the elastic Green function. 
The region 1 and 2 estimates are computed using the folded dynamical matrix in cubic supercells, and 
combined as in Eqn. (fTTl to produce the total deviation estimate, (a) Single random dynamical matrix shows 
an individual example of error estimation, (b) 100 different random dynamical matrices were computed, 
along with their associated LGF's. The average deviation over the ensemble average shows that we have 
an accurate computation of the error, even for the case of small supercells (2x2x2) with the long-range 
dynamical matrix (cutoff at 25ao). 

example, if the magnitude of the largest lattice vector i?niax is greater than 10/eias, the lattice Green 
function can be computed for lattice vectors |^| < lO/gias, and the elastic Green function used for the 
remainder, while introducing a total error of 1%. This choice can greatly speed the computation 
of the lattice Green function for large simulations by (1) limiting the ^-point grid size, and (2) 
restricting the set of points over which the full lattice Green function must be computed. 

The length scale /eias is also a fundamental length scale for quasi-continuum'*" and flexible 
boundary conditions methods^i^ where it determines the range at which the relaxation response 
using elastic finite-elements or the bulk continuum is accurate compared to atomistic response. 
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This length is not necessarily the same as the interaction force cutoff — it may be larger or smaller. 
For example, the region 2 estimate of deviation for an isotropic nearest-neighbor interaction gives 
^eias = V 1 / 6i?nn * OAR^n, which suggcsts transitioning from atomistic to finite-elements at twice 
the interaction cutoff produces errors on the order of 4% in position. On the other extreme, density- 
functional theory calculations in metals have shown surprisingly small /eias, considering the known 
long-range interactions in metallic systems.- It is the small value of Zgias that has allowed the ac- 
curate calculation of isolated dislocations using flexible boundary condition methods in density- 
functional theory. Knowledge of /gias is essential to constructing accurate computational cells that 
are large enough to produce accurate response, but do not waste computational resources treating 
interactions that can be replaced with elastic response. 

This paper presents an accurate computational algorithm for the lattice Green function from 
limited dynamical matrix information together with the elastic constants. In conjunction, an accu- 
rate error estimate using the limited dynamical matrix computable from a single supercell compu- 
tation allows measurement of the supercell-size effect. The error estimate produces a length scale 
4ias which characterizes the crossover from atomistic harmonic response to continuum elastic re- 
sponse. The algorithm for lattice Green function computation together with the determination of 
crossover length scale has already allowed the accurate computation of single extended disloca- 
tion defects using density-functional theory .^i^i^i"^^ The approach can also be utilized to implement 
flexible boundary condition methods for point defects, crack opening and tip propagation, surfaces 
and boundaries coupled with density-functional theory: providing chemically accurate interactions 
coupled with correct treatment of the long-range elastic response of extended defects. 
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APPENDIX A: ANGULAR INTEGRATION IN INVERSE FOURIER TRANSFORM 

The integration of the angular portion of the inverse Fourier transform in three dimensions can 
be performed analytically. Particular integrals referenced below can be found in Gradshteyn and 
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Ryzhik^^ 

Evaluation of the three dimensional integral 



JJ d'k e^^'^'-'^YUkX (Al) 



An 



begins by rotating the variable of integration to a new coordinate system given by p{k) such that R 
aligned along the /j^-axis. The spherical harmonic Yi,n{k) is written as a linear expansion of Yi^'ip) 

I 

YM = a^L^R)yim'{p), (A2) 

m'=-l 

where the coefficients will be determined later. With this expansion, Eqn. (lAlD becomes 

Val,(^) sinV'^^™^^" d<PpYj^,{ep,(l)p\ 



where {9p, (f>p) are the angular coordinates of p with R as the z-axis. In this coordinate system, 6p 
is the angle between p and R; hence, cos^p = k ■ R. The integral simplifies by (1) transforming 
u = cos 6p, and (2) noting that the (f>p integral is non-zero only for m' = 0. Thus, our integral 
reduces to 



a^^'l^iR) yfn V2Z + 1 ^ du P,{u)e''''^\ 



The integral of the Legendre polynomial are expressions 7.393.1 and 7.393.2 in Gradshteyn and 
Ryzhik^^, which has 

du Pi{u)e''" = 2{iiji{x), 



where ii{x) is the regular spherical Bessel function,— y/njlx Ji+i/2(x). 

Finally, a^^^iR) must be evaluated to produce the final expression. Eqn. iA2h can multiplied by 
Yi^(p) and integrated over 4n to give 

af!n(R) = Jj d^p YiMPWlip) 



An 

j2J 



= Jjd'k YMYl^ipik)). 

An 



Then, the n = component is 



a^^/R) = JJ d^k YUk) ^J^^P,(cos Op), 



An 
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where 6p is the angle between p(k) and the z-axis in p's coordinate system. Given our rotation, 
cos 6 = k ■ R. The addition theorem for spherical harmonics, 



m'=-l 



then gives 



m'=-l . 

4;r 



4;r 
2Z+ 1 



Combining the terms gives 



^ e'^^'^^^-^^yz^C^) = 4n(i)'ji(kR)YUR). 

An 



(A3) 



APPENDIX B: REGION 2 ERROR ESTIMATE 

The region 2 error estimate in Eqn. ([T6b contains two summations over ^ with _D(x) which can 
approximated using the / = spherical harmonic components of and G^'^ . The quartic term 
Yjx^^Wk^ appears in the computation of G^. 



Gdc 
00 



= G\h Y^x-h'm 



An 



G^{k). 



We can approximate G^{k) with its spherical average value, G^q/ ^^An. The integral ff^d^k (x-k)'^ 
4nx^/5, so 



J]x'Di£) 



120^[G^,S'g'^[G^]~\ 



The term requires more egregious approximations. Starting with the definition of Gqq, 



An 



both sides are inverted, and the inverse of the spherical average is approximated as the average of 
the inverse. 



L J? 



An 

1 An\ 
4^T2 



32 



Inverting again gives 

X 

Then, for the large R limit of Eqn. ([T6b . we have 



1 



6 V47r 



'00- 



1 

6^ 



10 G''^ 



00 



^00 



(Bl) 



The final approximation is to limit the summations to |^ < R; for that evaluation, we replace the 
true Gqq with Gqq(7?), evaluated using summations over |^ < R. 

Despite the approximations at use in Eqn. (IBll) . it is exact in an important limit: elastically 
isotropic crystals. Since the elastic Green function is isotropic, G^(k) = Gqq/ V47r, and the approx- 
imations in inverting Gqq are exact. Furthermore, for nearest-neighbor interactions, the left-hand 
side of Eqn. (|Bl|l is exactly R\J6R^, making /eias = Vl/6i?nn- 
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3nr for a review of phonons in ionic crystals. 
4^ This is the static bulk lattice Green function. The lattice Green function can be defined as a function of 
both time and space to model phonon propagation. In addition a defect lattice Green function can be 
computed using the Dyson equation. There, the starting point is the bulk lattice Green function. 
In fact, the sum of the lattice Green function over all lattice sites diverges; this requires that the sum of 
all forces on an infinite crystal body must vanish. 
4^ It should be noted that the elastic Green function possess a pole at the origin; hence, the deviation 
functional are technically infinite at ^ = ±R. This complication is avoided by choosing finite values for 
and e'^^ at the origin. This gives finite values for h,^^\G^^^W, +R) and A<2)[G^](^, +R). While this 
issue is ignored in the derivation to avoid complicating the notation further, it does not affect the results, 
as the elastic Green function is only evaluated at lattice sites in our derivation. 
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